Code
train_length <- floor(0.9 * length(inflation_monthly))
train_series <- inflation_monthly[1:train_length]
test_series <- inflation_monthly[(train_length + 1):length(inflation_monthly)]ARMA/ARIMA/SARIMA Models
Summary
Forecasting Inflation with ARMA/ARIMA/SARIMA Models
Inflation forecasting is of paramount importance for policymakers, investors, and consumers alike. Accurate predictions can aid in monetary policy decisions, investment strategies, and household budgeting. To forecast inflation using time series data, several models can be employed, notably ARMA, ARIMA, and SARIMA models.
Autoregressive (AR) and Moving Average (MA) Models:
AR: This model predicts future inflation based on its own past values. It assumes that the inflation rate is a linear function of its previous values.
MA: Contrary to AR, the MA model predicts inflation based on past white noise or error terms. This captures the shock effects observed in the past.
Autoregressive Moving Average (ARMA) Model:
- Combines both AR and MA components. Suitable for time series data that exhibit patterns not captured by AR or MA models alone.
Autoregressive Integrated Moving Average (ARIMA) Model:
An extension of the ARMA model that includes “integration” (I). This represents the number of differences needed to make the time series stationary.
Especially pertinent for inflation data, which may have trends or cycles that render the data non-stationary.
Seasonal Autoregressive Integrated Moving Average (SARIMA) Model:
- Extends ARIMA by accounting for seasonality, which can be pivotal for inflation data affected by seasonal factors (e.g., holiday-driven consumer spending or agricultural harvest cycles).
Checking Stationarity: For accurate forecasting, it’s vital that the inflation data is stationary, meaning its statistical properties like mean and variance remain constant over time.
Use the Autocorrelation Function (ACF) plot to test for stationarity. If significant correlations persist over several lags, the data may be non-stationary.
Differencing the data, often first or even second differences, can help in achieving stationarity by eliminating trends or cyclical patterns.
In our EDA section, an ACF plot was generated to examine stationarity in the inflation data. After ensuring stationarity, either naturally or through transformations, the ACF and Partial Autocorrelation Function (PACF) plots of the stationary data are pivotal in determining the optimal parameters for our ARIMA or SARIMA models.
Arima Modelling for Inflation
Splitting the data into Train and Test for Model Validation Process
Upon completing the cleaning and aggregation process of the U.S. Monthly Inflation Data, we are now poised to split this time series dataset into distinct training and testing subsets. This segregation is vital to ensure our forecasting model is both trained on a comprehensive data set and validated against a segment of data it has not previously encountered.
The data spans from 1913 to 2023 for every month. We split the data into a 70% train and 30% test.This allocation strategy will allow us to forecast and validate our model’s performance on inflation trends spanning over years.
ACF and PACF Plots for Monthly Inflation
Code
train_series %>%
ggtsdisplay(main="ACF and PACF Plots of Monthly Inflation")
ADF Test for Monthly Inflation
H0: The time series is non-stationary. In other words, it has some time-dependent structure and does not have constant variance over time.
v/s
H1: The time series is stationary. In other words, it has no time-dependent structure and does have constant variance over time.
Code
adf.test(train_series)Warning in adf.test(train_series): p-value smaller than printed p-value
Augmented Dickey-Fuller Test
data: train_series
Dickey-Fuller = -5.5685, Lag order = 10, p-value = 0.01
alternative hypothesis: stationaryBecause the p-value from the ADF test is less than α = 0.05, we reject the null hypothesis and conclude that the monthly inflation series is stationary.Although the ADF states that the original series is stationary, the ACF plots, which clearly indicate seasonality and trend, are more reliable than the ADF test. Therefore, it is safe to conclude that the series non-stationary as per the ACF section above.
Log Transformation of Monthly Inflation and its first and second order diferrencing
Code
lx = log(train_series + abs(min(train_series)) + 0.01) # since inflation has negative values
dlx = diff(lx)
ddlx = diff(dlx, 12)
x = train_series
# Plot original series
plot.ts(x, main="Original Series")
Code
# Plot log-transformed series
plot.ts(lx, main="Log-transformed Series")
Code
# Plot first difference of log-transformed series
plot.ts(dlx, main="First Difference of Log-transformed Series")
Code
# Plot second difference of first difference
plot.ts(ddlx, main="Second Difference of First Difference")
Code
par(mfrow=c(2,1))
monthplot(dlx); monthplot(ddlx)
Simply taking log of the number of monthly infaltion does not make it stationary. First-differencing the log number of monthly inflation does, however, make the series stationary and this series should be employed for building our time series model. Keep in mind that because first-differencing was enough to make the series stationary, we do not need to second-difference it, helping us avoid over differencing the number of monthly inflation
ADF Test for Log First Order Differencing Monthly Inflation
H0: The time series is non-stationary. In other words, it has some time-dependent structure and does not have constant variance over time.
v/s
H1: The time series is stationary. In other words, it has no time-dependent structure and does have constant variance over time.
Code
adf.test(dlx)Warning in adf.test(dlx): p-value smaller than printed p-value
Augmented Dickey-Fuller Test
data: dlx
Dickey-Fuller = -16.463, Lag order = 10, p-value = 0.01
alternative hypothesis: stationaryBecause the p-value from the ADF test is less than α = 0.05, we reject the null hypothesis and conclude that the log first-differenced monthly inflation series is stationary. Let us now check whether the ACF plots supports this hypothesis.
ACF and PACF Plots of Log First-Differenced Monthly Inflation
Code
dlx %>%
ggtsdisplay(main="ACF and PACF Plots of Log First-Differenced Monthly Inflation")
p values obtained from PACF are 0, 1, 2, 3, 4 q values obtained from ACF are: 0, 1 d (Difference): 1
Fitting ARIMA (P,D,Q)
Code
d=1
i=1
temp= data.frame()
ls=matrix(rep(NA,6*24),nrow=24) # roughly nrow = 3x4x2
# p=0,1,2,3,4
# q=0,1,2,3 (although we only found q=1 to be significant in ACF, we may want to compare a complex ARIMA model with greater "q" value compared to a simpler ARIMA model)
# d=1
{
for (p in 1:5)
{
for(q in 1:4)
for(d in 1)
{
if(p-1+d+q-1<=8)
{
model<- Arima(lx,order=c(p-1,d,q-1))
ls[i,]= c(p-1,d,q-1,model$aic,model$bic,model$aicc)
i=i+1
#print(i)
}
}
}
}
temp= as.data.frame(ls)
names(temp)= c("p","d","q","AIC","BIC","AICc")
temp <- temp[order(temp$BIC, decreasing = FALSE),]
knitr::kable(temp)| p | d | q | AIC | BIC | AICc | |
|---|---|---|---|---|---|---|
| 7 | 1 | 1 | 2 | -46.40664 | -26.0630351 | -46.37303 |
| 11 | 2 | 1 | 2 | -45.49580 | -20.0662959 | -45.44534 |
| 8 | 1 | 1 | 3 | -45.37218 | -19.9426758 | -45.32172 |
| 15 | 3 | 1 | 2 | -43.56630 | -13.0508907 | -43.49559 |
| 12 | 2 | 1 | 3 | -42.64466 | -12.1292475 | -42.57395 |
| 19 | 4 | 1 | 2 | -41.40796 | -5.8066500 | -41.31360 |
| 18 | 4 | 1 | 1 | -34.62703 | -4.1116249 | -34.55633 |
| 6 | 1 | 1 | 1 | -18.32448 | -3.0667737 | -18.30433 |
| 14 | 3 | 1 | 1 | -27.86882 | -2.4393084 | -27.81835 |
| 3 | 0 | 1 | 2 | -16.92151 | -1.6638090 | -16.90136 |
| 20 | 4 | 1 | 3 | -41.36933 | -0.6821168 | -41.24791 |
| 10 | 2 | 1 | 1 | -21.00188 | -0.6582723 | -20.96826 |
| 2 | 0 | 1 | 1 | -10.53575 | -0.3639459 | -10.52568 |
| 4 | 0 | 1 | 3 | -18.58800 | 1.7556010 | -18.55439 |
| 16 | 3 | 1 | 3 | -23.97866 | 11.6226453 | -23.88431 |
| 17 | 4 | 1 | 0 | 52.89847 | 78.3279760 | 52.94893 |
| 13 | 3 | 1 | 0 | 72.04419 | 92.3877970 | 72.07780 |
| 9 | 2 | 1 | 0 | 122.29070 | 137.5484074 | 122.31085 |
| 5 | 1 | 1 | 0 | 218.52444 | 228.6962410 | 218.53451 |
| 1 | 0 | 1 | 0 | 493.69542 | 498.7813227 | 493.69877 |
| 21 | NA | NA | NA | NA | NA | NA |
| 22 | NA | NA | NA | NA | NA | NA |
| 23 | NA | NA | NA | NA | NA | NA |
| 24 | NA | NA | NA | NA | NA | NA |
Best Model in terms of AIC:
Code
temp[which.min(temp$AIC),] p d q AIC BIC AICc
7 1 1 2 -46.40664 -26.06304 -46.37303Best Model in terms of AICc:
Code
temp[which.min(temp$AICc),] p d q AIC BIC AICc
7 1 1 2 -46.40664 -26.06304 -46.37303Best Model in terms of BIC:
Code
temp[which.min(temp$BIC),] p d q AIC BIC AICc
7 1 1 2 -46.40664 -26.06304 -46.37303Model summary and error metrics of ARIMA(2, 1, 1):
Code
fit <- Arima(lx, order=c(1,1,2)) # no drift included
summary(fit)Series: lx
ARIMA(1,1,2)
Coefficients:
ar1 ma1 ma2
0.8797 -1.6559 0.6576
s.e. 0.0462 0.0687 0.0664
sigma^2 = 0.05595: log likelihood = 27.2
AIC=-46.41 AICc=-46.37 BIC=-26.06
Training set error measures:
ME RMSE MAE MPE MAPE MASE
Training set 0.002588435 0.2361338 0.1155882 0.6293244 14.72852 0.8532608
ACF1
Training set 0.01477834The best model with the lowest BIC metric is the ARIMA(1,1,2).
AR(2): This indicates an autoregressive term of order 1. It means the model uses the two most recent values (lags) of the time series in the prediction equation. The AR(1) component suggests that there’s a relationship between an observation and the one preceding observations.
I(1): This represents the integrated component of order 1. It indicates that the time series has been differenced once to make it stationary. Differencing is the process of subtracting the current value from the previous value. In the context of inflation or economic data, this can be thought of as the change in value from one period to the next.
MA(2): This is a moving average term of order 2. The MA component accounts for the relationship between an observation and a residual error from an MA process applied to lagged observations. In simpler terms, it’s a way to model the error of the prediction as a linear combination of error terms from the recent past.
The ARIMA(1,1,2) model can be represented by the following equation:
\[ [(1 - 0.885B)(1 - B)X_t = (1 + 1.6633B - 0.6648B^2)Z_t] \]
Where:
- \(( B )\) is the backshift operator, which represents a lag of 1 period. \(( B^2 )\) would represent a lag of 2 periods, and so on.
- \(( X_t )\) is the time series value at time \(( t )\).
- \(( Z_t )\) is the white noise error term at time \(( t )\).
- The coefficients 0.8931, 1.6770, and -0.6786 are the estimated parameters for the AR(1), MA(1), and MA(2) terms, respectively, based on the provided output.
- The term \(( (1 - B) )\) represents the differencing of order 1.
Model Diagnostics of ARIMA(1,1,2)
Code
model_output <- capture.output(sarima(lx, 1,1,2))
Code
cat(model_output[71:102], model_output[length(model_output)], sep = "\n")
Coefficients:
ar1 ma1 ma2 constant
0.8922 -1.6761 0.6761 1e-04
s.e. 0.0242 0.0409 0.0408 1e-04
sigma^2 estimated as 0.05564: log likelihood = 27.85, aic = -45.69
$degrees_of_freedom
[1] 1191
$ttable
Estimate SE t.value p.value
ar1 0.8922 0.0242 36.9201 0.0000
ma1 -1.6761 0.0409 -40.9654 0.0000
ma2 0.6761 0.0408 16.5827 0.0000
constant 0.0001 0.0001 1.4324 0.1523
$AIC
[1] -0.03823631
$AICc
[1] -0.03820818
$BIC
[1] -0.01695639
NA
NA
NA
NA
NAStandardized Residuals: Essentially stating that if the errors are white noise. The model does look stationary as it captures all the signals and essentially captures the raw white noise.
ACF Of Residuals: Auto-correlation of the residuals. The only q value to inspect is 1.
Q-Q Plot: The series follows a normal distribution pretty closely as even the tails seem to be on the normal line.
p values of the Ljung-Box statistic: Ideally, we would like to fail to reject the null hypothesis. That is, we would like to see the p-value of the test be greater than 0.05 because this means the residuals for our time series model are independent, which is often an assumption we make when creating a model. Since the lag value less than 10 have a p-value greater than 0.05, the residuals have no autocorrelations.
Let’s see what model is outputted by auto.arima().
Model Output of Log Monthly Inflation with auto.arima()
Code
fit = auto.arima(lx, seasonal = FALSE)
summary(fit)Series: lx
ARIMA(2,1,1)
Coefficients:
ar1 ar2 ma1
0.1676 0.0859 -0.9026
s.e. 0.0442 0.0408 0.0317
sigma^2 = 0.05723: log likelihood = 14.5
AIC=-21 AICc=-20.97 BIC=-0.66
Training set error measures:
ME RMSE MAE MPE MAPE MASE
Training set 0.0003750843 0.2388288 0.1160352 1.793008 15.99104 0.8565603
ACF1
Training set -0.004293837From the above output, auto.arima() too outputted an ARIMA(2,1,1) model.Some points to keep in mind when using these functions is as follows:
The auto.arima() function in R uses a stepwise algorithm to search through the space of possible ARIMA models and select the one with the lowest AIC value. While this approach can be computationally efficient and provide a good starting point for model selection, it does not necessarily always find the best possible model for a given time series.
On the other hand, the Arima() function in R allows us to specify the exact order of the ARIMA model and can be used to fit more complex models, such as those with seasonality, exogenous variables, or other constraints. By specifying the exact order of the model, we have more control over the modeling process and can potentially obtain a better fit to the data.
In summary, the auto.arima() function can be a useful tool for quickly identifying a potentially good model, but it is not a substitute for careful model selection and customization seen when using the Arima() function.
The ARIMA(1,1,2) model outperforms the ARIMA(2,1,1) model for the lx series, as evidenced by its lower AIC, BIC, and log likelihood values, indicating a better balance between model fit and complexity. The ARIMA(1,1,2) model’s significant autoregressive term suggests the data has inherent sequential dependencies, which this model captures more effectively than the solely moving average-based ARIMA(2,1,1). While both models have similar error metrics, the simpler ARIMA(1,1,2) provides a more parsimonious and better-fitting representation of the underlying data dynamics.
Forecasting ARIMA(1,1,2) and ARIMA(2,1,1)
Code
arimaModel_1 <- arima(lx, order = c(1,1,2))
arimaModel_2 <- arima(lx, order = c(2,1,1))
forecast1=predict(arimaModel_1, length(test_series))
forecast2=predict(arimaModel_2, length(test_series))
# Convert the time series and forecast objects to data frames
ts_df <- data.frame(date = time(inflation_monthly), value = as.numeric(inflation_monthly))
train_df <- data.frame(date = time(inflation_monthly)[1:train_length], value = as.numeric(lx))
forecast1_df <- data.frame(date = time(inflation_monthly)[(train_length + 1):length(inflation_monthly)], value = forecast1$pred)
forecast2_df <- data.frame(date = time(inflation_monthly)[(train_length + 1):length(inflation_monthly)], value = forecast2$pred)
# Plot the time series and forecasts
ggplotly(ggplot() +
geom_line(data = train_df, aes(x = date, y = value,
color = "Actual Train Values"), linetype = "solid", alpha=0.6, show.legend = TRUE) +
geom_line(data = forecast1_df, aes(x = date, y = value,
color = "ARIMA(1,1,2) Forecast"), linetype = "solid", show.legend = TRUE) +
geom_line(data = forecast2_df, aes(x = date, y = value,
color = "ARIMA(2,1,1) Forecast"), linetype = "solid", show.legend = TRUE) +
geom_line(data = ts_df[(train_length + 1):length(inflation_monthly),], aes(x = date, y = log(value),
color = "Actual Forecast Values"), linetype = "solid", show.legend = TRUE) +
labs(x = "Date", y = "Log of Number of Monthly Inflation", title = "Forecasting ARIMA(1,1,2) and ARIMA(2,1,1)") +
theme_minimal() +
scale_color_manual(name = "Forecast",
values = c("ARIMA(1,1,2) Forecast" = "blue",
"ARIMA(2,1,1) Forecast" = "green",
"Actual Forecast Values" = "orange",
"Actual Train Values" = "black"),
labels = c("ARIMA(1,1,2) Forecast",
"ARIMA(2,1,1) Forecast",
"Actual Forecast Values",
"Actual Train Values")))From the above graph, we can note that the forecasted monthly inflations remains constant at around 1 for both models on the test set. This performance is not what was expected and, hence, it is possible that the models are not able to capture the underlying patterns in the data. This can be due to a variety of reasons, such as insufficient data and the models not being complex enough to capture the variation in the data. It is, however, pragmatic to check whether the sarima.for() function’s predictions may forecast differently. Let us find out below.
Forecasting Arima (1,1,2) using Sarima.for()
Code
log_monthly_inflation <- ts(lx, start = c(1913, 1), frequency = 12)
sarima.for(ts(train_df$value, start = c(1913, 1), frequency = 12), 24, p = 1, d = 1, q = 2, main = "Forecasting ARIMA(1,1,2) using sarima.for()") 
$pred
Jan Feb Mar Apr May Jun Jul Aug
2012
2013 1.266403 1.268101 1.269626 1.270997 1.272230 1.273340 1.274341 1.275243
2014 1.278622 1.279123 1.279581 1.279999 1.280382 1.280734 1.281058 1.281358
Sep Oct Nov Dec
2012 1.257422 1.260048 1.262401 1.264511
2013 1.276059 1.276796 1.277464 1.278071
2014
$se
Jan Feb Mar Apr May Jun Jul
2012
2013 0.2517885 0.2538952 0.2555667 0.2568957 0.2579543 0.2587988 0.2594735
2014 0.2614751 0.2616205 0.2617380 0.2618330 0.2619101 0.2619726 0.2620235
Aug Sep Oct Nov Dec
2012 0.2359764 0.2414628 0.2457539 0.2491269
2013 0.2600132 0.2604456 0.2607923 0.2610708 0.2612947
2014 0.2620650 Comparing Arima(1,1,2) with Benchmarks
Code
autoplot(log_monthly_inflation) +
autolayer(meanf(log_monthly_inflation, h=24),
series="Mean", PI=FALSE) +
autolayer(naive(log_monthly_inflation, h=24),
series="Naïve", PI=FALSE) +
autolayer(snaive(log_monthly_inflation, h=24),
series="SNaïve", PI=FALSE)+
autolayer(rwf(log_monthly_inflation, h=24, drift=TRUE),
series="Drift", PI=FALSE)+
autolayer(forecast(Arima(log_monthly_inflation, order=c(1,1,2)), 24),
series="ARIMA(1,1,2)",PI=FALSE) +
guides(colour=guide_legend(title="Forecast")) +
ylab("Log of Monthly inflation") + ggtitle("Forecasting ARIMA(1,1,2) and Benchmarks") + theme_minimal()
Arima(1,1,2) model metrics:
Code
fit <- Arima(log_monthly_inflation, order=c(1,1,2))
summary(fit)Series: log_monthly_inflation
ARIMA(1,1,2)
Coefficients:
ar1 ma1 ma2
0.8797 -1.6559 0.6576
s.e. 0.0462 0.0687 0.0664
sigma^2 = 0.05595: log likelihood = 27.2
AIC=-46.41 AICc=-46.37 BIC=-26.06
Training set error measures:
ME RMSE MAE MPE MAPE MASE
Training set 0.002588435 0.2361338 0.1155882 0.6293244 14.72852 0.7438047
ACF1
Training set 0.01477834Mean metrics:
Code
f1 <- meanf(log_monthly_inflation, h=24)
checkresiduals(f1)
Ljung-Box test
data: Residuals from Mean
Q* = 815.6, df = 24, p-value < 2.2e-16
Model df: 0. Total lags used: 24Code
accuracy(f1) ME RMSE MAE MPE MAPE MASE
Training set 2.290299e-16 0.2628595 0.1357713 3.834076 21.53009 0.8736818
ACF1
Training set 0.361126Snaive metrics:
Code
f2 <- snaive(log_monthly_inflation, h=24)
checkresiduals(f2)
Ljung-Box test
data: Residuals from Seasonal naive method
Q* = 579.05, df = 24, p-value < 2.2e-16
Model df: 0. Total lags used: 24Code
accuracy(f2) ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.0001674226 0.3444004 0.1554013 5.803652 24.34328 1 0.2192816Random Walk metrics:
Code
f3 <- rwf(log_monthly_inflation, h=24)
checkresiduals(f3)
Ljung-Box test
data: Residuals from Random walk
Q* = 335.5, df = 24, p-value < 2.2e-16
Model df: 0. Total lags used: 24Code
accuracy(f3) ME RMSE MAE MPE MAPE MASE
Training set 0.0001084129 0.2972421 0.1354665 3.373002 17.91581 0.8717203
ACF1
Training set -0.4551488The ARIMA(1,1,2) model seems to be the best performing among the four models considered. It has a MASE value less than 1, suggesting better performance than a naive model, and its residuals show low autocorrelation, indicating that the model has effectively captured the underlying patterns in the data. In contrast, the other models, especially the Mean and Seasonal Naive models, have not performed as effectively, with their residuals showing higher autocorrelation. The Random Walk model’s negative ACF1 suggests possible over-differencing. Therefore, for future forecasts and analyses on this dataset, the ARIMA(1,1,2) model would be the recommended choice. The performance of the ARIMA(1,1,2) model versus the other models can be attributed to the inherent nature of the models and the underlying patterns in the data:
Nature of ARIMA:
- ARIMA (AutoRegressive Integrated Moving Average) models combine autoregressive (AR) and moving average (MA) elements, allowing them to capture both autocorrelations and moving average patterns in the data. The (1,1,2) order suggests that the model uses one lag of the autoregressive term, one order of differencing, and two lags of the moving average term.
Simplicity of Other Models:
The Mean model uses a simple average of all values, which won’t capture any intricate patterns in a time series beyond its central tendency.
The Seasonal Naive model assumes a repetitive pattern after a fixed number of periods. If the data doesn’t have a consistent seasonality or has other complex trends, this model won’t capture them effectively.
The Random Walk model is based on the premise that future values are the same as the last observed value. This model is too simplistic for series with trends or seasonality, as it doesn’t anticipate any change in direction.
Dynamics of the Dataset:
- The errors (residuals) of a good model should be random without any discernable pattern. The ARIMA(1,1,2) model had the lowest autocorrelation in residuals among all the models, indicating that it has effectively captured most of the structure from the data.
Model Diagnostics:
- The metrics provided (AIC, BIC, MASE, ACF1) offer different perspectives on model performance. AIC and BIC values give information about the goodness of fit of the model, penalizing for complexity. MASE is a scale-free error metric that compares the given model to a naive benchmark, and ACF1 checks for any remaining autocorrelation in the residuals.
Monthly Inflation SARIMA Modelling
Code
# Visualize Seasonal Component
inflation_ts <- ts(inflation_monthly, frequency = 12) # frequency depends on your data, 12 for monthly
inf.s=decompose(inflation_ts)$seasonal
plot(inf.s, axes=FALSE,
main='Seasonal Component of Monthly Inflation in the US Over Time', xlab="Time", ylab = "Inflation", type='c')
Quarters = c("1","2","3","4")
points(inf.s, pch=Quarters, cex=1, font=4, col=1:4)
axis(1, 1:4); abline(v=1:4, lty=2, col=gray(.7))
axis(2); box()
From the above seasonal component graph of the number of monthly inflation, we notice there does exist some level of seasonality in the original series. The seasonal component graph illustrates the degree of seasonal variation in the monthly inflation. The magnitude of the seasonal variation is shown on the y-axis of the graph, and it indicates how much the monthly inflation deviates from the average value for each season. The graph shows a repeating pattern in the number of monthly inflation over time, with clear peaks in the first and second quarters and troughs in the third quarter. This pattern implies that the monthly inflation in the US might be influenced by the season of the year.
Seasonally Differenced Monthly Inflation
Code
# Seasonal differenced
attacks.diff=diff(inflation_ts,12)
plot(attacks.diff, axes=FALSE, main='Monthly Inflation (S. differenced)',type='c', ylab = "Inflation") #with type='c' I get dashed lines
Quarters = c("1","2","3","4")
points(attacks.diff, pch=Quarters, cex=1, font=4, col=1:4)
axis(1, 1:4); abline(v=1:4, lty=2, col=gray(.7))
axis(2); box()
ACF and PACF Plots of Seasonally Differenced Monthly Inflation
Code
inflation_ts %>%
diff(lag=12) %>%
ggtsdisplay(main="ACF and PACF Plots of First Seasonal Differenced Monthly Inflation")
After first ordinary differencing the original series (ACF?), we saw a lot of seasonal correlation left, suggesting that first order differencing did not help in transforming the raw data into a stationary series. This differenced series cannot be used for building a robust SARIMA model. Therefore, a seasonal differencing on the original monthly inflation was performed above and we can still notice some correlation left, but lesser compared to when the raw series was differenced with first order. Therefore, it could be that D=1 and d=0. Let’s keep this as one option and let’s proceed with performing both seasonal differencing and first-order differencing the raw monthly inflation series.
ACF and PACF Plots of Seasonally and First Ordered Monthly Inflation
Code
inflation_ts %>%
diff(lag=12) %>%
diff() %>%
ggtsdisplay(main="ACF and PACF Plots of Seasonally and First Order Differenced Monthly Inflation")
After both seasonal differencing and ordinary differencing together the raw data, the ACF and PACF plots seem to portray the least correlation than the individual differencing methods. Next, we shall difference and select the relevant p,d,q,P,D,Q values from the original monthly inflation series for our SARIMA model.
From the seasonal differencing and ordinary differencing (together) ACF and PACF plots, the following combinations for p,d,q,P,D,Q are:
q values obtained from ACF = 0,1,2,3,4 Q values obtained from ACF = 1 p values obtained from PACF = 0,1,2,3,4 P values obtained from PACF = 1,2 d (Difference) = 1 D (Seasonal Difference) = 1
Code
######################## Check for different combinations ########
#write a funtion
SARIMA.c=function(p1,p2,q1,q2,P1,P2,Q1,Q2,data){
#K=(p2+1)*(q2+1)*(P2+1)*(Q2+1)
temp=c()
d=1
D=1
s=12
i=1
temp= data.frame()
ls=matrix(rep(NA,9*29),nrow=29)
for (p in p1:p2)
{
for(q in q1:q2)
{
for(P in P1:P2)
{
for(Q in Q1:Q2)
{
if(p+d+q+P+D+Q<=9) # parsimonious principle
{
model<- Arima(data,order=c(p-1,d,q-1),seasonal=c(P-1,D,Q-1))
ls[i,]= c(p-1,d,q-1,P-1,D,Q-1,model$aic,model$bic,model$aicc)
i=i+1
#print(i)
}
}
}
}
}
temp= as.data.frame(ls)
names(temp)= c("p","d","q","P","D","Q","AIC","BIC","AICc")
temp
}
# q=0,1,2,3,4; Q=1 and PACF plot: p=0,1,2,3,4; P=1,2; D=1 and d=1
output=SARIMA.c(p1=1,p2=5,q1=1,q2=5,P1=1,P2=3,Q1=1,Q2=2,data=inflation_ts)
#output
knitr::kable(output)| p | d | q | P | D | Q | AIC | BIC | AICc |
|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 0 | 0 | 1 | 0 | 3480.647 | 3485.829 | 3480.650 |
| 0 | 1 | 0 | 0 | 1 | 1 | 2639.797 | 2650.161 | 2639.806 |
| 0 | 1 | 0 | 1 | 1 | 0 | 3076.969 | 3087.334 | 3076.978 |
| 0 | 1 | 0 | 1 | 1 | 1 | 2641.209 | 2656.757 | 2641.228 |
| 0 | 1 | 0 | 2 | 1 | 0 | 2973.552 | 2989.099 | 2973.570 |
| 0 | 1 | 0 | 2 | 1 | 1 | 2643.163 | 2663.893 | 2643.194 |
| 0 | 1 | 1 | 0 | 1 | 0 | 2988.168 | 2998.533 | 2988.178 |
| 0 | 1 | 1 | 0 | 1 | 1 | 2135.257 | 2150.804 | 2135.275 |
| 0 | 1 | 1 | 1 | 1 | 0 | 2584.681 | 2600.228 | 2584.699 |
| 0 | 1 | 1 | 1 | 1 | 1 | 2135.320 | 2156.049 | 2135.350 |
| 0 | 1 | 1 | 2 | 1 | 0 | 2454.577 | 2475.306 | 2454.608 |
| 0 | 1 | 2 | 0 | 1 | 0 | 2978.301 | 2993.849 | 2978.320 |
| 0 | 1 | 2 | 0 | 1 | 1 | 2124.433 | 2145.163 | 2124.464 |
| 0 | 1 | 2 | 1 | 1 | 0 | 2580.931 | 2601.660 | 2580.961 |
| 0 | 1 | 3 | 0 | 1 | 0 | 2939.600 | 2960.329 | 2939.630 |
| 1 | 1 | 0 | 0 | 1 | 0 | 3212.549 | 3222.913 | 3212.558 |
| 1 | 1 | 0 | 0 | 1 | 1 | 2388.030 | 2403.577 | 2388.048 |
| 1 | 1 | 0 | 1 | 1 | 0 | 2789.768 | 2805.315 | 2789.786 |
| 1 | 1 | 0 | 1 | 1 | 1 | 2389.776 | 2410.505 | 2389.806 |
| 1 | 1 | 0 | 2 | 1 | 0 | 2690.159 | 2710.888 | 2690.189 |
| 1 | 1 | 1 | 0 | 1 | 0 | 2941.689 | 2957.236 | 2941.708 |
| 1 | 1 | 1 | 0 | 1 | 1 | 2124.106 | 2144.835 | 2124.136 |
| 1 | 1 | 1 | 1 | 1 | 0 | 2580.136 | 2600.866 | 2580.167 |
| 1 | 1 | 2 | 0 | 1 | 0 | 2901.388 | 2922.118 | 2901.419 |
| 2 | 1 | 0 | 0 | 1 | 0 | 3125.210 | 3140.757 | 3125.228 |
| 2 | 1 | 0 | 0 | 1 | 1 | 2286.999 | 2307.728 | 2287.029 |
| 2 | 1 | 0 | 1 | 1 | 0 | 2718.776 | 2739.505 | 2718.806 |
| 2 | 1 | 1 | 0 | 1 | 0 | 2918.961 | 2939.691 | 2918.992 |
| 3 | 1 | 0 | 0 | 1 | 0 | 3013.639 | 3034.368 | 3013.669 |
Code
cat("\n Best Model in terms of AIC: \n")
Best Model in terms of AIC: Code
output[which.min(output$AIC),] p d q P D Q AIC BIC AICc
22 1 1 1 0 1 1 2124.106 2144.835 2124.136Code
cat("\n Best Model in terms of AICc: \n")
Best Model in terms of AICc: Code
output[which.min(output$AICc),] p d q P D Q AIC BIC AICc
22 1 1 1 0 1 1 2124.106 2144.835 2124.136Code
cat("\n Best Model in terms of BIC: \n")
Best Model in terms of BIC: Code
output[which.min(output$BIC),] p d q P D Q AIC BIC AICc
22 1 1 1 0 1 1 2124.106 2144.835 2124.136The best model with the lowest AIC, AICc and BIC metric is the SARIMA(1,1,1,0,1,1) model. The equation of the SARIMA(1,1,1,0,1,1) model is given by:
\[ \begin{align*} (1 - \Phi B^s)(1 - B^s)Y_t &= (1 - \theta B)(1 - B)X_t \\ \text{Where:} \\ Y_t &\text{ is the differenced series (after applying non-seasonal and seasonal differencing).} \\ X_t &\text{ is the original time series.} \\ B &\text{ is the backshift operator. For example, } B^k X_t = X_{t-k}. \\ \theta &\text{ is the parameter for the non-seasonal moving average component.} \\ \Phi &\text{ is the parameter for the seasonal autoregressive component.} \\ s &\text{ is the seasonal period.} \end{align*} \] ## Model Diagnostics of ARIMA(1,1,1)(0,1,1)
Code
model_output <- capture.output(sarima(inflation_ts, 1,1,1,0,1,1,12))
Standardized Residuals: Essentially stating if the errors are white noise. The model does look stationary as it captures all the signals and essentially captures the raw white noise.
ACF Of Residuals: However, looking at the ACF of the Residuals gives us a definitive answer to whether the model is stationary. Because some spikes are not within the significance limits, the model is not being able to capture all the signal in the data.
Q-Q Plot: The series weakly follows a normal distribution as the tails waver away significantly from the normal line.
p values of the Ljung-Box statistic: Ideally, we would like to fail to reject the null hypothesis. That is, we would like to see the p-value of the test be greater than 0.05 because this means the residuals for our time series model are independent, which is often an assumption we make when creating a model. Since all lag values greater than 5 have a p-value less than 0.05, residuals have remaining autocorrelations.
Forecast for the next 3 years using ARIMA(1,1,1)(0,1,1)
Code
fit <- Arima(inflation_ts, order=c(0,1,1), seasonal=c(0,1,1))
summary(fit)Series: inflation_ts
ARIMA(0,1,1)(0,1,1)[12]
Coefficients:
ma1 sma1
-0.7628 -0.9416
s.e. 0.0210 0.0093
sigma^2 = 0.2897: log likelihood = -1064.63
AIC=2135.26 AICc=2135.27 BIC=2150.8
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set -0.0001692118 0.5351804 0.3408826 NaN Inf 0.7206196 0.07360051Code
fit %>% forecast(h=36) %>% autoplot() #next 3 years
Comparing (1,1,1)(0,1,1) with Benchmarks
Code
cat("Best model metrics: \n")Best model metrics: Code
fit <- Arima(inflation_ts, order=c(1,1,1), seasonal=c(0,1,1))
autoplot(inflation_ts) +
autolayer(meanf(inflation_ts, h=36),
series="Mean", PI=FALSE) +
autolayer(naive(inflation_ts, h=36),
series="Naïve", PI=FALSE) +
autolayer(snaive(inflation_ts, h=36),
series="SNaïve", PI=FALSE)+
autolayer(rwf(inflation_ts, h=36, drift=TRUE),
series="Drift", PI=FALSE)+
autolayer(forecast(fit,36),
series="fit",PI=FALSE) +
guides(colour=guide_legend(title="Forecast"))
Code
cat("Best model metrics: \n")Best model metrics: Code
summary(fit)Series: inflation_ts
ARIMA(1,1,1)(0,1,1)[12]
Coefficients:
ar1 ma1 sma1
0.1367 -0.8271 -0.9450
s.e. 0.0375 0.0232 0.0092
sigma^2 = 0.2868: log likelihood = -1058.05
AIC=2124.11 AICc=2124.14 BIC=2144.84
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set -0.000610487 0.5323361 0.3387781 NaN Inf 0.7161705 0.001247957Code
cat("Snaive metrics: \n")Snaive metrics: Code
f2 <- snaive(inflation_ts, h=36)
accuracy(f2) ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.0008879703 0.7759247 0.473041 NaN Inf 1 0.3160833Seasonal cross validation of ARIMA(1,1,1)(0,1,1) and (2,1,1) using 1 Step Ahead Forecasts
Code
k <- 75 # minimum data length for fitting a model
n <- length(inflation_ts)
#n-k # rest of the observations
set.seed(123)
farima1 <- function(x, h){forecast(Arima(inflation_ts, order=c(1,1,1),seasonal=c(0,1,1)), h=h)}
e <- tsCV(inflation_ts, farima1, h=1)
MAE1 <-abs(mean(e,na.rm=TRUE))
cat("MAE for ARIMA(1,1,1)(0,1,1) is: ", MAE1)MAE for ARIMA(1,1,1)(0,1,1) is: 0.326548Code
RMSE1=sqrt(mean(e^2, na.rm=TRUE)) #one-step time series cross-validation
cat("\nRMSE for ARIMA(1,1,1)(0,1,1) is: ", RMSE1)
RMSE for ARIMA(1,1,1)(0,1,1) is: 0.716712Code
farima2 <- function(x, h){forecast(Arima(inflation_ts, order=c(2,1,1),seasonal=c(0,1,1)), h=h)}
e <- tsCV(inflation_ts, farima2, h=1)
MAE2 <-abs(mean(e,na.rm=TRUE))
cat("\nMAE for ARIMA(2,1,1)(0,1,1) is: ", MAE2)
MAE for ARIMA(2,1,1)(0,1,1) is: 0.32673Code
RMSE2=sqrt(mean(e^2, na.rm=TRUE)) #one-step time series cross-validation
cat("\nRMSE for ARIMA(2,1,1)(0,1,1) is: ", RMSE2)
RMSE for ARIMA(2,1,1)(0,1,1) is: 0.716795Both MAE and RMSE metrics agree that ARIMA(1,1,1)(0,1,1) is the best model by a slight margin. However, the BIC metric does not agree with this result as it outputted ARIMA(0,1,1)(0,1,1) as the model with lowest BIC. AIC and AICc metrics, however, do agree with the MAE and RMSE metrics generated from Seasonal Cross Validation using 1 step ahead forecasts. Let’s see whether this is the case when forecasting 12 steps ahead.
Seasonal Cross Validation of ARIMA(1,1,1)(0,1,1) and ARIMA(2,1,1)(0,1,1) Using 12 steps (Seasonal Period) Ahead Forecasts
Code
k <- 75 # minimum data length for fitting a model
n <- length(inflation_ts)
n-k # rest of the observations[1] 1254Code
set.seed(123)
farima1 <- function(x, h){forecast(Arima(inflation_ts, order=c(1,1,1),seasonal=c(0,1,1)), h=h)}
# Compute cross-validated errors for up to 12 steps ahead
e <- tsCV(inflation_ts, forecastfunction = farima1, h = 12)
mse1 <- colMeans(e^2, na.rm = TRUE)
farima2 <- function(x, h){forecast(Arima(inflation_ts, order=c(2,1,1),seasonal=c(0,1,1)), h=h)}
# Compute cross-validated errors for up to 12 steps ahead
e <- tsCV(inflation_ts, forecastfunction = farima2, h = 12)
# Compute the MSE values and remove missing values
mse2 <- colMeans(e^2, na.rm = TRUE)
# Plot the MSE values against the forecast horizon
data.frame(h = 1:12, MSE1 = mse1, MSE2 = mse2) %>%
ggplot() + geom_point(aes(y=MSE1,x= h)) + geom_point(aes(y=MSE2,x= h)) +
geom_line(aes(y=MSE1,x= h,colour="MSE for ARIMA(1,1,1)(0,1,1)")) +
geom_line(aes(y=MSE2,x= h,colour="MSE for ARIMA(2,1,1)(0,1,1)"))+
theme_minimal()
This plot gives cross-validation statistics up to horizon 12. The procedure for seasonal cross validation using 12 steps ahead is very similar to seasonal cross validation using 1 step ahead. We need to change the “h” parameter to the desired the number of time horizons we want to forecast for. The farima() function manually written by us helps us call our desired SARIMA model with the number of horizons. Then, farima() function is called inside the tsCV() function, which helps us store the cross-validated errors for up to 12 steps ahead. Then, because we get forecasts for each time horizon, we need to take the mean of the squared column using colMeans to obtain MSE.
Although we observed that the MSE and RMSE of ARIMA(1,1,1)(0,1,1) when forecasting 1 step ahead was lower than that of ARIMA(2,1,1)(0,1,1), from the above plot it can be seen that the cross-validated MSEs get lower or better as the number of forecasting steps increases. Both models’ MSE performance follow a very similar pattern, with ARIMA(1,1,1)(0,1,1), picked by lowest BIC, having a lower MSE across all forecasting steps, except for step 1. Therefore, ARIMA(2,1,1)(0,1,1) is the better SARIMA model!